First-order phase transition of fixed connectivity surfaces 
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We report a numerical evidence of the discontinuous transition of a tethered membrane model 
which is defined within a framework of the membrane elasticity of Helfrich. Two kinds of phantom 
tethered membrane models are studied via the canonical Monte Carlo simulation on triangulated 
fixed connectivity surfaces of spherical topology. A surface model is defined by the Gaussian term 
and the bending energy term, and the other, which is tensionless, is defined by the bending energy 
term and a hard wall potential. The bending energy is defined by using the normal vector at each 
vertex. Both of the models undergo the first-order phase transition characterized by a gap of the 
bending energy. The phase structure of the models depends on the choice of discrete bending energy. 

PACS numbers: 64.60.-i, 68.60.-p, 87.16.Dg 



I. INTRODUCTION 

Tethered membrane models 0, S IE are ordinar- 
ily defined by Hamiltonian that is a linear combination 
of discrete bending energy and surface tension energy 
[MIS' Hence, there may be a variety of statistical models 
of membranes, since discrete Hamiltonian can be chosen 
arbitrarily even within Helfrich or Polyakov-Kleinert pre- 
scription of membranes. As a consequence , it is natural 
to ask whether the phase structure [3, la, 111 III! LUj of the 
model depends on the Hamiltonian. 

However, little attention has been given to the depen- 
dence of the phase transition on the Hamiltonian of teth- 
ered surfaces both for models that have surface tension 
[3 [H Q 111 EJHJllJllJp EH and for tension- 
less models [2I IH X2A l25l l26l \2l\ . Almost all numerical 
studies done so far utilize the bending energy of the or- 
dinary form 1 — rij • n.j , where is the normal vector of 
the triangle i. 

One other discrete bending energy that has been uti- 
lized by Gompper et.al 28] is based on the discretization 
of the Laplacian in the dual lattice formulation of discrete 
mechanics by T.D.Lee [29] . Similar discrete bending en- 
ergy was adopted in Refs. [Ij, [13, y(J]- Both discrete 
bending energies give results compatible with the contin- 
uous phase transition of the model 0, H, 0, El EH • 

Recently, it was reported [3lJ that a tethered mem- 
brane model with the ordinary bending energy under- 
goes the discontinuous phase transition predicted in [33 , 
although the L-J potential is assumed to serve as the 
Gaussian term. Hence, we think it is worthwhile to show 
that the discontinuous phase transition can be seen in a 
tethered membrane model when the Hamiltonian is de- 



fined only by a discretization of Helfrich Hamiltonian. 

The purpose of this paper is to show numerical evi- 
dence that the phase structure of phantom tethered mod- 
els depends on the choice of the discrete bending energy. 
By using the normal vector at each vertex, we define a 
bending energy which is different from the ordinary one. 
We will study two kinds of models; one is a model that 
has the Gaussian term for surface tension and the other 
is a tensionless model that has no surface tension term 
but has a hard wall potential. It will be shown that both 
models undergo the first-order phase transition. 



II. MODEL AND MC TECHNIQUES 

Membrane models are ordinarily defined by the dis- 
crete Hamiltonian S = Si+feS^ with the bending rigidity 
6, where <Si is the surface tension energy and S2 is the 
bending energy respectively defined by 



Si — 2J [Xi - Xj 

(ij) 



s 2 = (i - n > 

(y) 



(1) 



E 



(ij) 



in Eq.JTJ is over all bonds (ij), and rij, n., are the 



unit normal vectors of the triangles sharing the bond (ij). 
Xi(£ R 3 ) in Si is the position of the vertex i. 

Other possible bending energy S2 can be obtained by 
using the normal vector of the vertex i such as 



n(i) 



Ni 

IN,; 



j(i) 



(2) 
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where Ej(i) denotes the summation over triangles j(i) 
linked to the vertex i. The vector n-jU) is the unit normal 
of the triangle j(i), and A^., n is the area of j(i). 
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The new discrete bending energy can be obtained by 
using the normal vector of Eq.|J2J. Thus, we have 



^ = EE[ 1 - n W-° 3 -w] 



(3) 



which is clearly different from that of Eq.Q. It should 
be noted that 5 2 (illdcf) = J2i j (1 — n (*) ' n (j)) defined 
only by using the normal vectors n(i) in Eq.(|2J) is not 
well defined. This ill-dcfinedness comes from the fact 
that there exists two surfaces locally different from each 
other that has the same value of 5*2 (illdef ) . Two normal 
vectors at the ends of a bond can be parallel for 

surfaces that are not smooth. 

We study two kinds of models in this paper. The first, 
which will be denoted by Model- 1, is a model defined by 

N 

Zi= \\dXieiq,[-{Si + bS 2 )}, (Model- 1) 

Sx = Y, (Pi - , 5 2 = E E ( x - n « • n m) > ( 4 ) 

(ij) • 

where the center of the surface is fixed to remove the 
translational zero mode. S2 is identical with J3J). 

The second model, which will be denoted by Model-2, 
is a tensionless model defined by 

N 

22 = / II dX * e M-(bS 2 + V)}, (Model-2) 
J i=i 

S a=EE( 1 - n W-n J -w). ( 5 ) 

where S% is identical with that of Model- 1 in Eq.Q, and 
V is the hard wall potential defined by 



V(\X i -X j \) = 



(0< \Xi-Xj] <ro), 
00 (otherwise) . 



(6) 



The value of r$ in the right hand side of © is fixed to — 
vL15. As a consequence we have (%2(Xi — Xj) 2 ) /N ~ 
3/2, which holds for Model-1 which contains the Gaus- 
sian term Si. It should be noted that Model-2 is consid- 
ered to be independent of the hidden length introduced 
by r . The MC results are independent of the value of 
tq. This was, in fact, precisely checked in Ref. |l5j . 

Figures ^ ( a ) an d Q (b) show the range of interactions 
described by S2 in © and the ordinary 52 in A 
difference between S2 in and S2 in 0} can be seen 
in the number of triangles whose normal vectors interact 
with the one of a given triangle, which is shaded in Figs. 
H(a) and^(b). The number of triangles for S2 in © is 
dependent on a given triangle and hence locally changes, 
while the number for S2 in is always 3. 

We use the canonical Metropolis Monte Carlo tech- 
nique. Spheres are triangulated by linking uniformly 
scattered points. The histograms of coordination number 
of surfaces are identical with those shown in Ref. 1151. 




FIG. 1: (a) Ranges of interaction between normal vectors of 
triangles for S2 in O), and (b) those for S2 in Q. The nor- 
mal vector of the shaded triangle interacts with those of the 
surrounding triangles in (a) and (b). Small spheres represent 
vertices. 



The position X of vertices is updated in MC by mov- 
ing the current position A" to a new position X' = X+5X, 
where ^A is chosen in a small sphere by using uniform 
random numbers. The radius Rq of the small sphere is 
fixed to Rq = do, where Iq is the mean bond length which 
is computed at every 250 MCS (Monte Carlo Sweeps), 
and a constant e is fixed at the beginning of the simula- 
tion to maintain 50 ~ 55% acceptance rate for Model-1 
and 55 ~ 65% for Model-2. The radius Rq becomes al- 
most constant, because is constant in the equilibrium 
configurations. 

We impose the lower bound 10~ 6 Aq on the area of tri- 
angles, where Aq is the mean area of triangles computed 
at every 250 MCS. As a consequence, the update of X 
are constrained so that the resulting area of triangles be- 
comes larger than 10 _6 ^4o- However, areas of almost all 
triangles are larger than 10~ 6 Aq in our MC without the 
lower bound, hence it seems that the areas are almost 
free from such constraint. No constraint is imposed on 
the bond length. 



III. RESULTS 
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FIG. 2: (a) Si/N B vs b of Model-1, and (b) Si/N B vs b of 
Model-2, where Nb is the total number of bond. iV = 1500. 

We first show Si/Nb of Model-1 and Model-2 respec- 
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tively in Figs. El (a) and|2(b), where Nb is the total 
number of bond. It should be noted that Si/Nb is the 
mean bond length squares Iq. Si/Nb in Fig. (a) of 
Model- 1 is completely compatible with the expected re- 
sult Si/N = 3/2, since iV B = 3iV-6(~ 37V) on the spher- 
ical surfaces. In fact, a typical sample in Fig. (a) is 
S' 1 /iV B = 0.50015±0.00012at 6 = 0.476. Moreover, Sx/N B 
in Fig. |21 (b) of Model-2 is also compatible with our ex- 
pectation Si/iV~3/2 as already stated in the paragraph 
below Eq. |JBJ|, although the Gaussian term Si is not 
included in the Hamiltonian of the Model-2. Thus, we 
confirmed that Iq is constant in the equilibrium configu- 
rations in both models. 

The specific heat Cs 2 is a fluctuation of the bending 
energy and is given by 



c S2 = w {(s 2 2 ) - (s 2 y 

Total number of MCS is about 0.< 



Cs2 
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(7) 
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FIG. 3: (a) Cs 2 vs b, and (b) Cf^ vs N in log-log scale, 
both (a) and (b) are obtained by Model-1 whose Hamiltonian 
is Si+bS 2 - (c) Cs 2 vs b, and (d) C™ 2 ax vs N in log-log scale, 
both (c) and (d) are obtained by Model-2 whose Hamiltonian 
is bS 2 + V. 

N = 340, TV = 600, 1.5 x 10 8 for N = 1000, 3 x 10 8 for 
N = 1500, and 2.2 x 10 8 for N = 2500 at the transition 
points b c (N) for Model-1. Number of MCS at b^b c (N) 
is relatively small. Total number of MCS for Model-2 is 
smaller than that for Model-1, since the speed of conver- 
gence of Model-2 is relatively faster than Model-1. 

Figure |3 (a) shows Cs 2 vs b of Model-1. The peak 
values C!? ax of Model-1 are plotted in Fig. 0(b) against 
N in log-log scale. Figures (c) and (d) are results 
obtained by Model-2. The number of molecules are N = 
340, TV = 600, TV =1000, iV=1500 for Model-2. 



The slope of the straight lines in Fig. 09(b) and Fig. |3| 

(d) represents the critical exponent a defined by 



rfmax 
°5 2 



N° 



(8) 



The largest three data in each figure are included in the 
fit, and we have 



a x = 0.798(66) 
a 2 = 0.914(166) 



(Model-1) 
(Model-2) 



(9) 



The value ax — 0.798(66) is smaller than 1 and hence 
implies that the order of the phase transition of Model-1 
is of second-order. However, as we will see next, the order 
of the phase transition of Model-1 is considered to be of 
first order. While ax < 1, the value a 2 =0.914(166) almost 
equals to 1 and hence suggests that Model-2 undergoes 
the first-order phase transition. 
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FIG. 4: (a) Variation of S2/NB against the number of MCS, 
and (b) the histogram h(S2 / 'N b) , obtained by Model-1. The 
results obtained by Model-2 are shown in (c) and (d). 

To clarify the order of the transition of Model-1, we 
plot in Fig. 01 (a) the variation of S 2 /Nb against the 
number of MCS. The series S 2 shown in Fig. (a) was 
obtained at every 5 x 10 4 MCS at the transition point b — 
b c (N) on the surface of size TV = 2500. The corresponding 
histogram Ii(S 2 /Nb) is drawn in Fig. E|(b). Figures 0] 
(c) and 0](d) are the results obtained by Model-2 of size 
TV = 1500. 

We clearly see in Fig. |U(a) that there are two distinct 
states which represent a discontinuous phase transition in 
Model-1. The histogram in Fig. 01(b) shows more clearly 
the existence of the two states separated by a gap of S2 
in Model-1. It is also easy to understand from Figs. 0](c) 
and (d) that Model-2 undergoes the first-order phase 
transition characterized by a gap of S 2 . 
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FIG. 5: (a) Variation of X 2 against the number of MCS, and 
(b) the histogram h(X 2 ), obtained by Model-1. The results 
obtained by Model-2 are shown in (c) and (d). 



The mean square size X 2 defined by 



X 2 



(10) 



is plotted in Fig. Efa) against the number of MCS of 
Model-1. The corresponding histogram h(X 2 ) is drawn 
in Fig. GJb). Figures (c) and (d) are the results 
obtained by Model-2. We see two different sizes at the 
transition point in each model, hence consider that the 
phase transitions in both models are characterized also by 
discontinuity of X 2 . The reason why we use X 2 obtained 
at b = 0.443 in Figs. El (c),(d) is that the double peaks 
in the histogram of X 2 at b = 0.443 is clearer than at 
6 = 0.442 where the histogram of S^/Nb plotted in Figs. 
El(c),(d) were obtained. 

The Hausdorff dimension [3^, |3 IHl is defined by 

X 2 ~N 2 ' H . (11) 

The gap of X 2 at the transition point implies that H 
discontinuously changes at that point. 

We plot in Fig. El (a) X 2 vs b of Model-1. The mean 
square size X 2 at b = b c (N) are plotted against N in 
log-log scale in Fig. EJb). The straight line denoted by 
smooth is obtained by fitting X 2 , each of which is the 
larger X 2 in the double peaks shown in Fig. El (b). An- 
other straight line denoted by crumpled is obtained by 
fitting the smaller X 2 in the peaks. Errors of X 2 were 
not included in the least squares fitting, since the fitting 
was done by using only the peak values of X 2 in the his- 
togram shown in Fig. 0(b). FiguresEl(c) andEl(d) show 
the results of Model-2. 
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FIG. 6: (a) X 2 vs b of Model-1, (b) X 2 vs N at the transition 
point b = b c (N). The results obtained by Model-2 are shown 
in (c) and (d). 



From the slope of the straight lines in Figs. El (b) and 
El (d), we have 



fli(T) = 2.13(17), Hi(i) = 3.66(107) 
fla(T) = 2.16(30), H 2 (i) = 7.84(977) 



(Model- 1), 
(Model-2). 



(121 



H(1)(H(l)) is considered as the Hausdorff dimension in 
the smooth (crumpled) phase at b>b c (b<b c ) just above 
(below) b c in each model. The reason of the large errors 
both in Hi (I) and in -Ha (J.) seems come from the fact 
that there is a few data points of X 2 included in the 
fitting. 

We understand from the straight lines in Figs. El (b) 
and El (d) that the phase transition of Model-2 is rela- 
tively stronger than that of Model-1, although both of 
the transition are the first order. The gap of H at b = b c 
of Model-2 is relatively larger than that of Model-1; this 
difference of H can be visible in the slope of the straight 
lines in Figs. E](b) andE|(d). 

There is no difference between the surfaces in the 
smooth phase of the models in this paper and those of 
[l5j . While the surfaces in the disordered (or crumpled) 
phase of the models in this paper are more crumpled 
than those in . The Hausdorff dimension at b > b c of 
the models and those of [15J are comparable, although 
the order of the transition of the models in this paper is 
different from that in [l^; the models in ^jjj have the 
continuous phase transition. 

We note also that both -Hi(T) and /^(T) are compat- 
ible with ( or slightly smaller than) the Flory prediction 
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H = 2.5, and they are almost compatible with an analyti- 
cal result H = 2.39(23) which corresponds to the scaling 
exponent z/ = 0.84±0.04 in ^ where v = 2/H. The val- 
ues -Hi(T) and -Had) m Eq. I|12l) imply that the surfaces 
are relatively swollen and smooth in the smooth phase at 
6>6 C in both models. 




FIG. 7: Snapshots of Model-1 surfaces obtained at (a) b = 
0.475(crumpled phase), (b) b = 0.478(smooth phase), and the 
sections of the surfaces in (a) and (b) are shown in (c) and 
(d) respectively. iV = 2500. 

In order to see the surfaces, we show snapshots of size 
N = 2500 of Model-1 in Figs. □ (a) and (b); one of 
which is obtained in the crumpled phase at 6 = 0.475 and 
the other in the smooth phase at 6 = 0.478. The sections 
of them are shown in Figs. (c), (d). The surface swells 



in the smooth phase as expected. We also find that the 
surface in Fig. (b) is smooth only at long range scales 
and rough at short scales. This is compatible with that 
seen in the model with the ordinary bending energy |15| . 
The surfaces of Model-2 are almost the same as those in 
Fig. El 



IV. SUMMARY AND CONCLUSIONS 

We have shown that the continuous phase transi- 
tion seen in the ordinary tethered membrane models is 
strengthened in two kinds of tethered membrane mod- 
els, whose bending energy is defined by using the normal 
vectors at the vertices. One of the models is defined by 
the Hamiltonian Si + bS2, and the other is a tensionless 
model defined by 6S2 + V, where V is a hard wall po- 
tential. It was shown by extensive MC simulations that 
both of the models undergo the first-order phase transi- 
tion which is characterized by a gap of S2. The size of 
spherical surfaces and the Hausdorff dimension discon- 
tinuously change at the phase transition in both models. 

The definition of the Hamiltonian remains in the 
framework of membrane elasticity of Helfrich. The bend- 
ing energy in Eq. © utilized in this paper appears to 
induce a non-nearest neighbor interaction between nor- 
mal vectors of the surface. In fact, the range of the inter- 
action is a bit larger than that of the ordinary bending 
energy as depicted in Fig. However, the bending en- 
ergy in Eq. @ is written by the normal vectors of Eq. 
(0) and the normal vectors of the neighboring triangles, 
and hence it is defined only by local geometric quantities 
of the surface just like the ordinary bending energy in 
Eq. 0. 
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